% Test some code in terms of simulated MLE
clear
clc
% person, period, drugid, choice, copay, prev


% ignore serial correlations
% focus on class 2

% sketch
% fminsearch @ (t) SMLE(t, data, random_draws)
% SMLE function: simulate using random draws scaled by parameters
% compute a total likelihood of seeing data as a function of parameters
% Compute a Hessian of the function at optimal parameters



% personid quarter drugid chosen cost_sharing incumbent control_extra treat_extra
% starts in June 1996 (1)
% can go to May 2006 (40)
A = importdata('general_entry_choice_data.txt');

data = A(:,1:6);
id = A(:,1);
quarter = A(:,2);
drugid = A(:,3);
chosen = A(:,4);
cost_sharing = A(:,5);
incumbent = A(:,6);

% % filtering
filter = (quarter <= 40) .* (drugid < 4) .* (1-isnan(cost_sharing)); % for now, filter based on brand, pre-2006
id(filter == 0, :) = [];
quarter(filter == 0, :) = [];
drugid(filter == 0, :) = [];
chosen(filter == 0, :) = [];
cost_sharing(filter == 0, :) = [];
incumbent(filter == 0, :) = [];
data(filter == 0, :) = [];


% extra parameters: number of periods, etc.
num_drugs = max(drugid);
num_periods = max(quarter);
num_people = max(id);

% setup work common to each loop
[~,~,id] = unique([id quarter], 'rows'); % unique person/quarter combos => creates choice id
outside = 1 - accumarray(id, chosen); % chose the outside option


% draw some normal random numbers
rng(17)
S = 100;
random_draws = normrnd(0,1,num_people*num_drugs, S);
random_draws2 = normrnd(0,1,num_people*num_drugs, S);
random_draws3 = normrnd(0,1,num_people*num_drugs, S);
random_draws4 = normrnd(0,1,num_people*num_drugs, S);

% parameters delta for each drug x period, alpha, gamma, sigma2 (of idiosyncratic)
% start with values that roughly reflect previous estimates
delta_init = [-1.0783;-0.3111;-1.6126];%[-0.82;-0.05;-1.35];
sigma_init = [0.0908;-0.9129;0.1589;-0.1346; -0.1; 1; 0.1; 1]; %[0.1;1;-0.1;1];
%params_init = repmat(delta_init, num_periods, 1); % -1*ones(num_drugs * num_periods, 1);
params_init = [delta_init; 0.0258; 5.2565; sigma_init];
params = params_init;

myopts = optimset('TolFun', 1e-10, 'MaxFunEvals', 12000, 'MaxIter',10000, 'display','iter'); % off, iter, final, notify

% test
SMLE_v2(params_init,data,random_draws,random_draws2,random_draws3,...
    random_draws4,num_drugs,num_periods,id,outside,S)

[params_opt, mle] = fminsearch (@(t)...
    -SMLE_v2(t,data,random_draws,random_draws2,random_draws3,...
    random_draws4,num_drugs,num_periods,id,outside,S), params_init, myopts);

save('smle_mixture4.mat');

delta = params_opt(1:3);

alpha = params_opt(4);
gamma = params_opt(5);

delta
[alpha, gamma]
params_opt(6:13) % mixture parameters



%% Compute standard errors
fun = @(t) -SMLE_v2(t,data,random_draws,random_draws2,random_draws3,random_draws4,num_drugs, num_periods, id, outside, S);


% hessian of a function at the final point
[H_small, H, H_big, G] = own_hessian_v3(fun, params_opt);

n=length(params_opt);
for i=1:n
    for j=i+1:n
        H(i,j) = H(j,i);
        H_small(i,j) = H_small(j,i);
        H_big(i,j) = H_big(j,i);
    end
end

variance_mat = inv(H_small);
standard_errors = diag(variance_mat.^(.5));

params_opt(4:13)
standard_errors(4:13)

save('smle_mixture4.mat');
